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Classical Taylor-Aris dispersion theory is extended to describe the transport of 
suspensions of self-propelled dipolar cells in a tubular flow. General expressions for 
the mean drift and effective diffusivity are determined exactly in terms of axial 
moments, and compared with an approximation a la Taylor. As in the Taylor-Aris 
case, the skewness of a finite distribution of biased swimming cells vanishes at long 
times. The general expressions can be applied to particular models of swimming 
microorganisms, and thus be used to predict swimming drift and diffusion in tubular 
bioreactors, and to elucidate competing unbounded swimming drift and diffusion 
descriptions. Here, specific examples are presented for gyrotactic swimming algae. 
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1. Introduction 

Suspensions of swimming microorganisms, such as algae and bacteria, behave dif- 
ferently to molecular fluids. Many microorganisms exhibit taxes, directed motion 
relative to external or local cues. For example, various algae (e.g. Chlamydomonas 
and Dunaliella sp.) swim upwards on average in the dark (gravitaxis) due either 
to a centre-of-mass offset from the centre-of-buoyancy (Kessler 1986), sedimenta- 
tion and anterior-posterior asymmetry in body/fiagella (Roberts 2006), or active 
mechanisms (Hader et al. 2005). This can result in aggregations of cells at upper 
boundaries and, if the cells are more dense than the medium in which they swim, 
overturning instabilities, termed bioconvection (Wager 1911; Piatt 1961). Further- 
more, a balance between gravitational and viscous torques can bias cells to swim 
towards downwelling regions, whence their added mass amplifies the downwelling. 
This is known as a gyrotactic instability and does not require an upper bound- 
ary. Of particular relevance here, Kessler (1986) observed that for a suspension 
of gyrotactic Ghlamydomonas nivalis in a vertically aligned tube, cells became 
sharply focused at the centre for downwelling flow and scattered towards the edges 
when the flow was upwelling. Additionally, phototrophic algae are often phototactic 
(they swim towards weak light and away from bright light), which can modify the 
instability mechanisms above, and bacteria may exhibit chemotaxis (e.g. up oxygen 
gradients). In shallow containers, the above taxes can result in very distinct bio- 
convection patterns, with characteristic lengthscales of millimeters to centimetres 
in just tens of seconds (Bees & HiU 1997). See Pedley & Kessler (1992) and HiU & 
Pedley (2005) for reviews. In deep cultures, one may observe long thin plumes of 
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Figure 1. Bioconvection plumes in suspensions of C. nivalis in a) a culture 
flask (5.3 cm wide), and b) a long vertical tube of inner diameter 1.1 cm. 
Concentration f» 1.5 x l(f cm~^. 



cells (Figure [1} that have a clear impact on the transmittance of light through the 
culture, of some relevance to photosynthetic algae (a "Cheese plant effect" ) . 

Recently there has been renewed interest in utilizing microorganisms for fuel 
production. For green algae there are two main approaches: hydrogen production 
by sulphur deficient cells (Melis & Happe 2002) and biomass generation for biodiesel 
production (Chisti 2007). To reach economical viability, both methods require the 
sustainable culture of cells, extensively and under carefully controlled conditions. 
Culture systems typically consist of arrays of tubes (vertical, horizontal or helical) 
and aim to maximize light whilst maintaining linear separation of cell stage and 
medium victuals. In algal bioreactors, suspensions of algae are typically pumped 
and may be bubbled or tansported turbulently to enhance nutrient/gas mixing 
and reduce variance in light exposure. These processes, which treat a suspension 
of microorganisms like a chemical fluid, are energetically costly. Instead, efflcient 
bioreactor designs might hope to harness the activity of the swimming microorgan- 
isms directly in laminar flows. However, it is unclear how a) the mean cell drift and 
b) the effective axial swimming dispersion of cells are affected by various flow flelds 
in the aforementioned tube arrangements. 

In a series of papers Taylor (1953, 1954a, 1954b) described how it is possible to 
approximate the effective axial diffusion of a solute in a fluid flowing through a tube. 
Molecular diffusion and advection by shear each play a distinguished role, such that 
the effective diffusivity is given by D„i + U'^a? / A%Dm, where Dm is the molecular 
diffusivity, a is the radius of the tube and U is the mean flow speed. Subsequently, 
Aris (1955) formalized the approach by solving the moment equations, extending 
ubiquitously the domain of physical relevance of Taylor's result. The methods have 
been extended by many authors (e.g. Horn & Kipp 1971, partitioning reactions 
between phases; Brenner 1980, dispersion in periodic porous media). The value of 
the Taylor- Aris approach can be measured by the wealth of practical applications 
(see Alizadeh et al. 1980). Until now, the approach has not been extended to sus- 
pensions of biased swimming microorganisms in a tubular flow. As we shall see, 
it is possible to derive general expressions with few assumptions. However, these 
expressions depend upon constitutive equations for the mean behaviour of the cells. 
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We shall adopt the standard continuum approach to modelling bioconvective 
phenomena, although our main result is independent of the details of these de- 
scriptions. Recent models of dilute, gyrotactic bioconvection (Childress et al. 1975; 
Pedley & Kessler 1990, 1992) assume that the fluid flow is governed by the Navier- 
Stokes equations with a negative buoyancy term to represent the effect of the cells 
on the fluid (Boussinesq approximation) such that 

= -Vpe + Apg + V • S, (1.1) 

where u(x, t) is the velocity of the suspension, Pe(x, t) is the excess pressure, I](x, t) 
is the stress tensor, g is the acceleration due to gravity, n(x, t) is the cell concen- 
tration, A/9 is the difference between the cell and fluid density, p, and v is the mean 
volume of a cell. The cell Reynolds number is small (e.g. ~ 10~^ for C. nivalis). 
Furthermore, the suspension is assumed incompressible such that V • u = 0. Pedley 
& Kessler (1990) extended the standard Newtonian description to include Batche- 
lor stresses, stress associated with rotary particle diffusion, and swimming induced 
stresslets. The first two were found to be qualitatively and quantitatively insignif- 
icant, and the third only plays a role in concentrated regions of the suspension. 
Thus in a dilute limit one may write V • S = /iV^u, where fi is the fluid viscosity. 
We shall employ this approximation in explicit examples, but the main result does 
not require it. Typically, over the course of a bioconvection experiment the total 
number of cells is conserved, so that one may write 

^ = -V-[n(u + Ve)-D.Vn], (1.2) 

where Vc(x) is the mean cell swimming velocity and D(x) is the cell swimming 
diffusion tensor, both of which need to be determined. At rigid boundaries, Q, we 
require a no-slip condition, u = on ^, as well as zero cell flux normal to Q (in 
direction n), such that n • (n (u + Vc) — D • Vn) = on ^. 

To model gyrotaxis, Pedley & Kessler (1987) employed a deterministic balance 
of gravitational and viscous torques on a spheroidal cell, of eccentricity ao; to 
determine the cell orientation p: 

p = ^ [k - (k • p)p] + ^HAp + aoP • E • (I - pp). (1.3) 

Here, B is the gyrotactic reorientation time-scale of a cell affected by external (grav- 
itational) torques subject to resisting viscous torques, given by B = iJ,a±/2hpg, 
where h is the centre-of-mass offset relative to the centre-of-buoyancy and a± is 
the dimensionless resistance coefficient for rotation about an axis perpendicular to 
p. n and E are the local vorticity vector and rate-of-strain tensor, respectively. 
These authors then wrote Vc = V^p, where Vs is the mean swimming speed, and, 
as for earlier models, assumed a constant isotropic diffusion. Pedley & Kessler 
(1990) advanced this description by postulating that the probability density func- 
tion, f{p,t), for orientation p satisfies a Fokker- Planck equation, with drift due to 
the various torques and a rotational diffusivity analogous to rotational Brownian 
motion (Frankel & Brenner 1991), thus taking account of biological variation of 
swimming stroke. Experimental data on cell tracking (Hill & Hader 1997) has pro- 
vided values for the deterministic and diffusive parameters. Prom /(p), the mean 
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swimming direction, q, is easily calculated, yielding Vc = Vsq, but the cell swim- 
ming diffusion tensor is not and requires approximation. Pedley & Kessler (1990) 
suggested that D « Vg^Tvar(p), where r is a direction correlation time, estimated 
from experimental data, and found asymptotic solutions for small flow gradients. 
Bees et al. (1998) extended these solutions for all flow gradients by expansion in 
spherical harmonics (employed in Bees & Hill 1998, 1999). However, the ad hoc 
nature of the diffusion approximation was cause for concern. This motivated Hill & 
Bees (2002) and Manela & Frankcl (2003) to develop generalized Taylor dispersion 
theory (Frankel & Brenner 1991), taking account of both the orientation and posi- 
tion of cells swimming in a linear flow, to derive the leading order, long time, spatial 
diffusion tensor. The techniques were subsequently employed by Bearon (2003) for 
dispersion of chemotactic bacteria in a shear flow. 

There are significant qualitative differences between the three treatments de- 
scribed above as vorticity is varied. In particular, as vorticity, uj, is increased the 
Fokker-Planck and the generalized Taylor dispersion approaches provide eigenvalues 
of the diffusion tensor that tend towards non-zero and zero limits, respectively. This 
is due to the fundamental difference between the orientation only versus trajectory 
based descriptions. Such qualitative differences in behaviour need to be tested with 
laboratory experiments. One approach is to track individual microorganisms in the 
very dilute limit (Hill & Hader 1997; Vladimirov et al. 2004) but for a precisely 
prescribed shear flow (e.g. Durham et al. 2009). However, such a scheme would 
likely be laborious and may not easily yield significant results for large shear rates. 
A macroscopic approach would be much preferred. In general, the coupling between 
cell and fiuid is bidirectional; the flow is driven by the presence of the cells, which 
determines the swimming directions of the cells. Controlling the flow in the manner 
described by Taylor may thus be advantageous. There are, however, some obstacles 
to be overcome. In particular, a local distribution of cells will drive secondary flows 
and lead to an effective axial diffusivity that depends on the axial location. The 
answer is to create a flow that is independent of the presence of the cells. This can 
be achieved by creating a long axisymmetric plume of swimming cells and dying a 
small blob of cells within the plume (Figure 2). In this way, we partially uncouple 
the drift-diffusive dynamics of the dyed cells from the bulk flow-cell problem. 

In the next section, we shall describe the geometry and scaling of the problem 
and introduce the method of moments. In section[3l steady-state solutions of plume 
concentration and flow in a tube subject to a pressure gradient are calculated. In 
section 01 the long-term drift and effective diffusion of a blob of cells in a plume 
in a tube of circular cross-section are formulated in general terms. The skewness of 
the distribution is also determined. For general comprehension and comparison, an 
argument in a vein similar to that given by Taylor (1953) is presented in section[5j 
The full theoretical results are then summarized in section[6]before explicit example 
calculations are given. Conclusions are presented in section [T] 

2. Flow in a straight tube 

This analysis is applicable to the case where the flow is independent of the axial 
direction. Thus consider the diffusion of dyed cells within a long plume. 

We follow the notation of Aris (1955) and consider a tube with characteristic 
scale a with axis parallel to the vertical a;-axis (pointing in the downwards direction; 
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Figure 2. Tube arrangement and coordinates. The magenta (light grey) region 
on the left represents dyed cells for which drift and diffusion is to be calculated. 



see Figured]). The interior of the tube is denoted by S, its cross-sectional area by s 
and its perimeter by F. We consider flows, u, generated by a pressure gradient and 
the added mass of the algae such that 

u{xh) = u{xH)ex = U[l + x{^H)]ex, (2.1) 

where U is the mean flow speed and x is the flow speed relative to the mean and 
is assumed to be only a function of the cross-sectional coordinates xh- Clearly, a 
no-slip boundary condition provides x ^1 on F. 

Let the cell swimming diffusion tensor be of the form D'^T), where D"^ is its 
characteristic scale, and the mean cell swimming velocity be Fsq(x//) (see Bees et 
al. 1998, where q = (p)), where Vs is the mean swimming speed. As x is independent 
of the axial direction then so are D and q. This fact permits a treatment using the 
method of moments in a similar vein to that described in Aris (1955). 

The cell conservation equation (jl.2p can thus be written 

^nt = V • (D • Vn) - (1 + x) - • (nq) , (2.2) 

where we use subscripts for partial differentiation where it is clear. It is conducive to 
translate to a reference frame travelling with the mean flow, and non-dimensionalize, 
such that X = [x — Ut)/a, = Xfl-/a and i ~ D'^t/a?. Equation (|2.2p becomes 

= V • (D • Vn) - Pexn, - /3V • (nq) , (2.3) 

where 

Pe=^, and = ^ ^ j , (2-4) 

and the hats are dropped for notational clarity. Here, Pe is a Peclet number, which 
is a ratio of the rate of advection by the flow to the rate of swimming diffusion, and 
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/3 is a 'swimming' Peclet number, a ratio of the rates of transport by swimming 
to swimming diffusion (or tube radius to swimming correlation length, where r is 
the direction correlation time, as typically = V^t\ e.g. Pedley & Kessler 1990). 
No-flow and no-flux boundary conditions shall be applied to the solution, such that 

u = and n • (D • Vn - /3qn) = 0, on T, (2.5) 

respectively, where n is normal to T. 

The pih moment with respect to the axial direction through xjj is defined as 

/ + 00 
xPn{x,XH ,t)dx, (2-6) 
-oo 

provided it exists and is finite (i.e. xPn{x,XH,t) — >■ as x — >■ ±oo). The cross- 
sectional average (denoted by an overbar) of this moment is written 

mp{t)—c^——J J CpdS. (2-7) 

Henceforth, consider axisymmetric fiows in a tube of circular cross-section with 
radius a oriented parallel to the vertical a:-axis (pointing downwards). Here x — 
x(r), D — D(r) and q = q(r) (such that q has no component in the eg direction). 

In cylindrical coordinates, by multiplying by x^ and integrating over the length 
of an infinite pipe, equation (12. 3p becomes 

Cp,t = ^ [?■(£)'''' Cp,r - (3q''cp - pl?™cp_i)]^ - pD™Cp_i,r 

+piPex + Pqncp-i +p(p- l)i?"%-2, (2.8) 

with 

D^'^'cp^r - Pq^'cp - pD^'^'cp^i =0 on r = 1. (2.9) 
Averaging over the cross-section (applying no-fiux boundary conditions 12.91) yields 

mp,t ^ -pD"Cp_i,^+p(Pex + /3g^)cp_i + p(p - l)D^^Cp_2. (2.10) 

Before deriving results for drift and diffusion in section 21 we shall solve the 
steady, coupled, cell conservation and hydrodynamic problem. 



3. Steady problem: flow and cell concentration 

Kessler (1986) demonstrated theoretically and experimentally that plume solutions 
exist in vertically aligned tubes. He found that the plumes are generally stable when 
a pressure gradient is applied such that the fiow is downwards. However, varicose 
instabilities may arise when no pressure gradient is applied. Here, we aim to avoid 
such instabilities and thus in the ensuing analysis implicitly refer to parameter 
regimes where plume solutions are stable. 

In later sections, in order to compute the dispersion of a blob of cells within a 
plume, we require knowledge of X: the fluid velocity relative to the mean. Hence, 
when x(^) represents the steady fluid velocity induced by a pressure gradient and 
the presence of a swimming cell distribution that is independent of x, 

= V-(D-Vn*)-^V-(qn*) (3.1) 
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where n* now represents all cells in the plume, and not just those dyed cells for 
which we shall calculate dispersion. As n* = = g'' at r = 0, this implies that 

D'-'-n;= pq'-n*. (3.2) 

Hence, given D^''{r) and q^{r) we have 

n = n(0) exp f /? r 4M^^ds] , (3.3) 



where fi is non-dimensional cell concentration (scaled with the average concentra- 
tion, n*). Note, for a spherical cell (ao = 0), and D^^ are functions of vorticity 
only, which must be in the eg direction: w = VAu = —Xri''')eg = iveg. 

In cylindrical polars the steady flow equation in the dilute limit becomes 

^ r ''^^''^r = - ah, (3.4) 

subject to the boundary conditions Xr(0) = and x(l) = —1. Here, the non- 
dimensional pressure gradient is Px = Pxo/UpL, and 



a^vgApN 



(3.5) 



measures the magnitude of the effect that the cells have on the flow, g is the 
acceleration due to gravity acting in the positive x-direction. 

Contrary to intuition, px and a are not free parameters but are linked to the 
mean flow speed, U, introduced in equation (j2.ip . Together they are determined by 
the boundary conditions on x and the requirement that x = 0? the flow deviation 
relative to the mean is order one. For Poiseuille flow, where a = 0, it is well known 
that Px = —8, such that x — 1 — 2r^. 

Substituting (13. 3p for n, equation (j3.4D can be rewritten as 

- {rXr)r -p. = -a m exp (p f 4^,^^) ■ (3-6) 



For spherical cells (ao — 0), taking logs and differentiating provides 



Note, differentiating removes the dependence on a; to fully specify the constants of 
integration, substitution back into Equation (13. 6p will be required. In general, equa- 
tion (j3.7p can be solved for uj and, thus, x and n (with application of the boundary 
conditions). Later, we shall consider the simple case 7(w) « Alj, for constant and 
negative A, so here we derive expressions for x in this limit. Equation p.7p becomes 

r'^uj" + (r - Aujr'^) lu' - {1 + rAuj)) uj = Pxr^Auj. (3.8) 

r = is a singular point and so consider u) = J2m=o brnf"^^^ , where constant Q is to 
be determined. Substituting into the nonlinear equation and examining coefficients 
reveals Q = I, for finite solutions at r = 0. Furthermore, the recurrence relation 



A 

bt = - 



Pxbt-2 + J2li=0 brnbt-rn-2{m + 2) 



tit + 2) 



(3.9) 
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is forthcoming. We require that w is odd and, therefore, hi = 0, Vi odd. Hence, the 
first few coefficients are given by 

b2 = ^ [Px + 26o] , 64 = [p^ + 26o] [Px + 660] , 

^6 = + + 66o]fe + 8feo] + ^[Px + 26o]'. (3.10) 

Furthermore, application of the boundary conditions yields 

x^-i + E^^d-^"^')- (3-11) 

m=0 

Applying the condition x = admits the result 

Finally, substitution of x into Equation p.6p is required to find a in terms of 
&,„, m = 0, 2, 4, and p^. Equation p. 61) can be written 

(^Xr)r - = -raexp [-Ax{r)], (3.13) 

where ci = aT%{0) exp (Ax(0)) (evaluated with the normalization condition fi — 
1, giving n(0)e'^'^*^°-' = 1/2 e~"^'^'^''V c?r). Hence, substituting x into p.l3p and 
comparing coefficients at leading order in r we find that 

5 = [2feo+P.]cxp|-A(^l- (3.14) 

Higher orders in r provide a check for the previously computed bm, m = 2,4,6, .... 
Therefore, given 6m, to = 0, 2, 4, and Px, then a can be computed from (j3.14p . 

If bo = 60(5) is required, and in the particular case that A is small (i.e. the cells 
are weakly affected by the flow; e.g. B is small), such that we can neglect quadratic 
terms in A and higher, then we can expand the transcendental equation to give 



-2p. + 25{l + A(l-E™=i;^)} 



bo ^ —■ (3.15) 

4 + aA 

Three examples are presented below, with the profiles plotted in Figure [3] 

(I) One of the simplest cases is for a = (i.e. the presence of the cells does 
not affect the flow). In this case, we compute bo = —px/'2 and bm — 0, 
TO = 2,4,6,..., such that x = ^1 + Pxi'r'^ — l)/4, which is Poiseuille flow. 
Equation p.l2p gives &o = 4 = —px/2, as would be expected. 

(II) With 5 7^ and A < small (i.e. a broad plume), but a zero pressure 
gradient, px — 0, then 62 = AbQ/A, bm = O(A^), to — 4,6,8,.... and equa- 
tions (irn)) and (I3J5)) provide bo = 6(-l ± ^/l + 8A/3)/2A + 0{A^) = 
<3 [1 + (1 - 5/4) A] + 0(^2). Thus X = -1 + boil - r^)/2 + Abl{l - r^)/16 + 



Article in press: Proceedings of the Royal Society A 



Dispersion of microorganisms in a tube 



9 




-1 -0.5 0.5 1 

Figure 3. Flow profiles for three cases (see text; extended to cover the width of 
the tube): a = 0, px ~ —8 (solid); a ^ 0, px = 0, A — —j (two dashed curves); 
a ^ 0, Px = —6, A — — J (dotted and dot-dashed curves). Selected broad cell 
distributions, n/ri(0)e'''"\ are also plotted (inset; see text). 

0{A^). Two solutions are possible: a simple positive flow (mode 1), and one 
with upwelling towards the edge of the tube (mode 2). For zero pressure gra- 
dient, a closed form, mode 1 solution is known. Kessler (1986) noted that 
h = n(0)/(l + Cin(O)r^)^ is a solution, for constant Ci. Applying the con- 
dition n = 1 gives Ci — n(0) — 1. Substituting this solution back into the 
governing equation reveals that n{0) is determined by the constrained pa- 
rameter a, as should be the case, in the same way that the mean velocity is 
linked to the pressure gradient in Poiseuille flow. This closed-form profile is 
approached by the above mode 1 profile with truncated sums (not shown). 

(Ill) The case A < (small; a broad plume) and p^; 7^ is also of interest, and has 
not previously been investigated. HA — —j and px — —6, then we calculate 
62 = -boibo - 3)/2\ 64 = boibo - 3)(&o - l)/2^ and feg = -6o(&o - 3)(&o - 
l)(46o - 3)/(2i33) - blibo - 3)2/(2^^3). From equation (|3A2)) we compute the 
two solutions 60 ~ 4.179 and 21.931. Again, the mode 2 solution corresponds 
to a flow with upwelling near the edge of the tube. The corresponding a 
can be evaluated from equation p.l4p . Hence, for the mode 1 solution, x ~ 
-l+4.179(l-r2)/2-0.308(l-r4)/4 + 0.0612(l-r6)/6-0.0107(l-r8)/8 + .... 

4. Dispersion in a tube of circular cross-section 

In this section, we place no restrictions on the cell shape (e.g. spheroidal cells) and 
form of x{r), q(r) and D(r), and find general expressions for the drift and effective 
diffusion of a blob of dyed cells within an existing plume. 

(a) Cell conservation and drift 

For p = equation (12.101) gives mo^t = 0, so that mp is a constant (i.e. number 
of cells is conserved). We fix mo = 1 (and remember that these cells represent dyed 
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cells diffusing within a plume of other cells) . Equation (|2.8p with p = implies that 

co,t = - [r{D^^CQ,r - /3g"co)], , (4.1) 
r 

with boundary condition D^'^co^r — f^q^co = on r = 1. The solution takes the form 

oo 

co = i?°(r) + 5]i?°(r)T°(i), (4.2) 

n=l 

where T° = exp (— 7,^jt) and satisfies 

r7?-i?0" + {D"- + rD^ i??/ + {-^q^ - r^q^' + j^r) = 0, (4.3) 
subject to the initial conditions. The solution for i?g, such that i?o(= cq) = 1, is 



Note also that ^ = 0, n 7^ 0. Putting p = 1 in (P?TU1) gives 

00 

TOi,t = Ao + X]cxp(-7f.i)A„, where A„ = -D™i?0' + (Pcx + /?g^) i?o . (4.5) 

ri=l 

In particular, it is clear that 



lim mi.f = Ao = -D^^R^' + (Pex + ^g^) i^n- (4.6) 

This means that the mean of the blob of dyed cells will move at a speed of Aq 
relative to the mean flow. Hence, 

mi(t) = Aot + ^ ^ (1 - exp (-7'*)) A„, (4.7) 

n=l 

where we have used mio = 0. The first term of Aq in (|4.6p is associated with a 
diffusive fiux, the second with advection of the cells heterogeneously distributed 
near the axis of the tube, and the third to swimming in the vertical direction 
relative to the fluid motion. At long times we expect 

00 ^ 

mioo(i) = Aoi + ^^, (4.8) 

ri=l 

(6) Effective diffusion 
With p = 1 , equation (j2.8l) implies that 

ci.t - - Wc^.r ' Pq^'ci - i?™co)], = -D-^c^^r + (Pe X + /3g")co, (4.9) 
r 

with boundary condition 

n'^'^ci,^ - /Jg'^ci - D^co = on r = 1. (4.10) 
The solution of this equation can be constructed in three parts. 
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1. Particular integral from i?o(r) in cq. It satisfies 



-ci, + - [r {D^^cl, - Pq^c\ - D^^Rl)] ^ = i?"i?o'-(Pex + m K (4.11) 

2. Particular integral from the rest of the terms R^{r) exp (— 7^0? 7^ 0, in cq. 
1 



- (r)'e-^"* - (Pox + /39") i?° (r)e-''"'. 



(4.12) 



It is quite clear that solutions to (|4.12l) are of the form S'„(r) exp (— 7^i), 
where Sn satisfy no-flux boundary conditions and are found by solving 



7^^„(r) + - [r (i?-5„(r)' - Pq^S{r) - D'''Rl{r))]^ 



D^-Rlir)' ~{Vex + mK{r)■ 



{A.r^) 



As we are interested in long-time behaviour, we do not solve for S'„(r), but 
later will require its cross-sectional average. This can be found by averaging 
both sides of (|4.13p and using the boundary conditions (|4.10l) . to give 



>-'n — 9 

7„ 



A. 



-D™i?o' + (Pex + /3g^)i?o = ^. (4.14) 



3. Complementary function. Solutions of (|4.9p without terms in cq that satisfy 
(j4.10p are of the form Al^R^{r)e~"'"* , where A^^ are constants. 

For item 1, to calculate c\, we rewrite the equation as 

cl,r- [r {D^^cl, - Pq^c\ - = r [-D^^R^ + (Pex + m 4) Mr), 

(4.15) 

RecaUing that Rl satisfies D'''' R^^ - Pq'' Rl = 0, let c\{r,t) = [Mt -t- /(r)] where 
M is a constant and /(r) is a function of r. Then, (j4.15p becomes 



[rU'D'^Rl - U'-Rl)]' = -Xo + MR°r, 
an equation independent of t. Hence, integrating once provides 



r(/'i^-i?° - i5™i?°) = -lA*{r) + ^Mm*ir) 



(4.16) 



(4.17) 



where 



A*{r) = 2 r Xo{s)ds = 2 r s + (Pex + /3g") <) ds, (4.18) 



mo{r) = 2 / sR'^{s)ds, 



(4.19) 

Ai5(l) = Ao and mQ{l) — 1. Applying the no- flux boundary condition (I4.10p to 
(|4.17p yields M = Aq. Integrating (|4.17p again provides 

c}(r, = i?[](r) (Aot + /(r)) = [Aq* + J(r) - $(r)] R°{r), (4.20) 
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where 



J(r) = 



, — — )-4-ds and <i>(r) = - 

D^'-{s) ^ ' 2 



sD^^{s)Rl[s) 

Hence, the complete solution ci — c\+ c\~\- c\ is given by 



ds. (4.21) 



ci = [Aot + J(r) - $(r)] i?[J + ^ 5„( 



re 



re 



(4.22) 



n=l 



n=0 



A\ are chosen to fit the initial data cw{r). In particular, the value of Aq is fixed 
by the initial condition mio = cTo = 0. With i?o = 0, n ^ 0, this implies 



- 

^0 — 



n=l 



, where J" = [$(r) - J(r)]i?o. 



(4.23) 



Thus the axial mean eventually is distributed across the tube as 



R°o{r) + [Jir)--i>{r)+^]R",ir) 



= {Al + Aot + f{r))R°{r). 



(4.24) 



After averaging across the cross-section we obtain cioo = mi<x> = AQt — J2n=i ^n{r), 
which can be compared with the earlier equation for the long-time limit of mi (j4.8p . 
and allows the identification 5„(r) = — A„/7^ consistent with equation (|4.14p . 

Putting p = 2, substituting the long-time solutions for ci and cq in (I2.10p . and 
using definition (|4.5I) for Aq, gives 



2D" [( J - $ + J-)i?E]]' + 2(Pex -I- I3q^) {J - ^ + J")i?[] 



-2An 



71=1 



(4.25) 



If De is the effective axial diffusion then one may define = limt^^oo 



where V is the variance {V = ^ Jg{x — xYndsdx). Then 

Dc = lim \^Am2-m\) 
t^oo 2 dt 



(4.26) 



= -D^^ [(J - $ + J')R°n]' + (Pex + Z??^) ( J - 3' + -7')i?[] + i'^^i?;] 



(c) Third moment and approach to normality 
With p = 3 equation (fTTO)) becomes 



mg^f -3i:»™C2,r + 3(Pex + /3(?^) C2 6D^^ci, 
where C2 is a solution to (|2.8p with p = 2, 

C2,t = i [r(i?'-''c2,. - /3g'-C2 - 2i^"ci)]^ - 2Z?'--ci, 
r 

+2(Pex + /3g")ci+2i^^^co, 



(4.27) 



(4.28) 
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subject to D^'^C2.r — Pq^C2 — 2D''^ci = on r = 1. The long time solution to 
equation (j4.28p has the form 

C2oo(r, t) = [2D,t + Alt' + 2AlAot + B^] i?E!(0+2(Aoi+^J)/(r)i?°(r)+.9(r)i?°(r), 

(4.29) 

where Bq is a constant determined by the initial distribution of dyed cells and 
function g(r) can be established after some algebra. Substituting (|4.24p and (|4.29p 
into KTi} gives 



m3,t 



Ao [2Det + Klt^ + 2AlK^t + Bl] + 2(Aot + Al){D^ - FK^) + H{r) (4.30) 



where H{r) ^ (gR^)' + (Pex + Pq'') gR^ + 2L>^^/i?[j and we have used the 

definitions (j4.6p and (j4.26p for Aq and respectively. Rearrangement and integra- 
tion thus provides 



ms - 3mim2 + 2m? = 3 H{r) - 2D^fRl - Aog{r)R^ 



t + const, 



(4.31) 



yielding the absolute skewness, y^, of the concentration distribution, with 



(to3 — 3miTO2 + 2mf ) 



3\2 9 



H{r) - 2DJRI - Kg[r)Rl 



(^200 - mf^)3 

(4.32) 

Hence, the skewness of the distribution decays to zero as i^^/^ as in classical Taylor- 
Aris dispersion; at long times we expect a Gaussian profile for the algal blob aver- 
aged across the cross-section. 



5. Mean drift and effective diffusion 'a la Taylor' 

It is instructive to re-derive approximations to equations (14. 6p and (j4.26p from (|2.3p 
using an approach similar to that of Taylor (1953). We use Taylor's approximations 
without a rigorous attempt to defend them. We begin by assuming that the cell 
concentration can be written as a superposition of the cross-sectionally averaged 
concentration, n — ri{x,t){= ^ J JgudS), given that it is well-defined, and a term 
Sn = Sn(x,r,t) for the radial variation, such that 

n{x,r,t) ~ n{x,t) + Sn{x^r,t). (5-1) 

Substituting (jS.ip into (|2.3p we find 

nt + 5nt = i {r[D''''Snr - f3q''{n + Sn) + ^^(rTj, + Sn^:)]}^ 

-l-D^^n™ - (Pex + I3q''){n^^ + Sn^) + D''''{n^^ + 5n^^) (5.2) 

subject to D^'^6nr — /3q^{n + 5n) + D''^{nx + Sn^) = on ?- = 1. Then, taking the 
cross-sectional average of both sides of (15. 2p gives 

nt = D™(5n^^ - P'^n^ - (Pex + I3q'')5n^ + D^n^x + D^^Sn.^^, (5.3) 



Article in press: Proceedings of the Royal Society A 



14 



M. A. Bees & O. A. Croze 



where we have used 5n = Q = x ^-nd the boundary condition. The aim is to express 
5n as a function of n to write (j5.3p in the form of an advection-difFusion equation 
for n. First, subtract (j5.3p from (15. 2p to obtain 



5nt = i{r[i?'^'-<5n,-/3<7'-(n + <5n) + i?'-^(n, + ^n,)]},, + /3^n, (5.4) 



-(Pex + Pq')na^ - [(Pex + Pq"')5n^ - (Pcx + Pq'')6n.^] + 



XX • 



Next, with Taylor (1953, 1954b), we make the following assumptions: (1) axial 
contributions to diffusion are negligible with respect to the radial ones and axial 
advection (V^n ^ V^n and Ux)', (2) concentration gradients in the axial direction 
are independent of radial position {Sux ~ 0; ); (3) transients decay rapidly; 

and (4), for simplicity, radial concentration fluctuations about the mean are small, 
5n <^ n. Note that the last assumption is not necessary, and is made only for 
illustrative convenience. Equation (|5.4[) then reduces to 

i {r[D"Snr - (3q^n + = [(Pex + /3g") - (5.5) 

subject to D^^Srir — fiq^'n + D^^rix = 0, on r = 1. Equation (|5.5p thus gives 

Sn = SRqU — {J ~ (f) + a)nx, (5-6) 
where SR'^ = j3 -^dr, and J — jy^dr, as in equation (|4.2ip . Furthermore, 
1 f 2/;a(Pex(a)+/3g-(a)~/3^)da 



(r) = - / ^ \. .\ ds, (5.7) 



and a = (0 — J) is a constant obtained by imposing Sn = 0. Using (j5.6p . (j5.3p reads 

nt + AqUx = DeUxx, (5.8) 
where we neglect terms of order Uxxx, consistent with previous approximations, and 



Ao = + (Pex + /3g=^)^i?§ + /3g^ (5.9) 



L»e = -D-'^{J - (j) + ay + {Pex + f3q'^)iJ - (l^ + a) + D'^^. (5.10) 

The above equations are limiting forms of (|4.6I) and (|4.26p . To see this, expand 
i?g ~ (1 + SRq), where SRq = (3 D^r^]^) ds <^ 1 (implying a broad distribution 
across the tube) . Substituting into (|4.6I) and (|4.26l) and neglecting terms of order 
{6Rq)^, leads to the above expressions. As earlier, there is a drift of cells relative to 
the flow due to swimming, diffusion and cell weighted average of the flow. 

6. Examples of dispersion 

(a) Summary of drift and effective 



To recap our main results, the drift, Aq, and effective axial diffusivity, D^, of 
a dyed blob of algae within an axisymmetric algal plume in a tube of circular 
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cross-section are given by 

Ao = + (Pex + /3g-)i?g, (6.1) 

D, = [(J - + (PCX + I3q- - Ao) (J - <P)R"o + (6.2) 

where Pe and (5 are Peclet numbers feauation l2.4p . 

To evaluate the above expressions we require the flow field relative to the mean, 
x(r), and constitutive equations for the mean cell swimming direction, q(r), and 
swimming diffusion tensor, D(r). Expressions for x(r) are obtained in section |3l 
and q(r), D(r) are available from solutions to deterministic or statistical models of 
gyrotaxis (Pedley & Kessler 1987, 1990; Bees et al. 1998; Hill & Bees 2002; Manela 
& Frankel 2003). 

The base distribution of cells, Roir), is defined by equation (|4.4p . Furthermore, 
the functions J{r) and <J>(r) are computed from equations (j4.2ip . and require the 
functions AQ^r) and mo(r) defined by equations (|4.18p and (|4.19p . respectively. 

(b) The limit to classical Taylor- Aris dispersion 

A useful check on the results is to reduce them to the original 'non-swimming' 
form of Taylor (1953) and Aris (1955). The original molecular solutes were assumed 
to diffuse isotropically (with no biased motion) and have no influence on the flow. 
Hence, put D'^^ = 1 = D^^, D^^ = (thus J{r) = 0) and qi = 0. For a circular pipe, 
Poiseuille flow provides x{i^) = 1 ^ 2r^. Thus, ^{r) = Pe^ /J ^ (/p ax{cr)da) ds — 

Pe — , so that x$ = Pe/48 and Rq = 1. Then the effective transport coeffi- 
cients (|6.ip and (|6.2p reduce to 

Ao and De = 1 + Pe^-?-, (6.3) 

48 

the classical Taylor- Aris result. In this same limit, equation (j4.24p for the the centre 
of mass of the solute distribution at long times reduces to cioo = —^{r) + $ + 
J2 ^=ii^ n/ln) SO that TOioo = = J2n=ii^ri / Jn) ^ whcrc $ = Pe/12 and A„ = 
Pex-R^. Thus cioo = mioo + Pe (l/12 - r'^/A + r*/8), consistent with Taylor- Aris. 

(c) Poiseuille flow limit for weak (rju ^ Ij and strong gyrotaxis (rjuj ^1^ 

As a second example, consider a simple Poiseuille flow not affected by the pres- 
ence of the cells. Then, x = 1 ~ 2r^ and uj = — Xr = 4r. We consider the limits of 
weak and strong gyrotaxis quantified by the ratio of the timescale for reorientation 
by the flow, = {U/a)~^, and the characteristic time-scale for reorientation of 
a cell by gravity against viscous resistance, B = 2mgh ' gyrotactic reorientation 
time. The ratio rj = BQ, is called the non-dimensional gyrotaxis parameter. Here, 
we consider two limits for which analytic solutions are known for the Fokker-Planck 
equation governing the probability distribution for the cell orientation p of spher- 
ical cells, due to Pedley & Kessler (1992) and Bees et al. (1998). Using definitions 
for the J and K constants from these papers, if 77 ^ 1 then — —Ki + 0{rf'uP'), 
q'- = -Jir?w-f 0(7^3^3)^ j^TT ^ Ki/X + Oirfuj^), D''^ = -riuj{J2~ JiKi) + 0{Tfuj^) 
and Z)^^ = K2 + 0{jf'uP'). At the other extreme, for 77 ^ 1 we have the asymptotic 
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solution = O {'q'^uj-^) , q'' = -l?r^uj-^ + O {t]~^uj-^) , D""" ^ i + O (ry-^^-^^ , 
L»" = O (rj-^uj-^) and D=''= = ^ + O {rj-^uj-'^). 

Substituting uj = Ar and omitting higher orders for clarity obtains, for <C 1, 

q'' = -4Ji77r, = -Ki, D"-^ = Girjr, D''-' = i^i/A, D''-' = X2, (6.4) 

where Gi — — 4( J2 — Jii^i), and, for -q^ 1, 

q"-^-- — , 9^ = = 13™, !)'■'■ = 1 =1?^^. (6.5) 
6 ?7r 3 

(i) Drift and effective diffusivity, rj <ti I 

In this limit the cells are less affected by the flow and prone to swim upwards. 
Using dm and defining = Ki/{2Ji\l3ri), equations (gjHl), (33) and 
provide 

A*(r) = mS(r) [2G177 + Pe(l - 2rl) - K^/S] + 2r^R° [Per^ - dv] , (6.6) 

which satisfies mQ{l) — 1, as required. Hence, in the limit rj ^ I the drift, Aq, is 

Ao = A5(l) = 2G177 (1 - i?Sj(l)) + Pe [1 - 2rl (l - R"o{1))] - (6.8) 

highlighting the contributions of swimming diffusion, advection and upswimming. 
In a similar manner, the expression (j6.2l) for the effective diffusivity becomes 

D, = -2Gnjao + [2Pcrl-2aiG2]Ii{l)-2Pch{l) + K2, (6.9) 

where oq = (J(l) - $(l))ai, ai = i?g(l), /,(r) = 2 /J's*(J(s) - $(s))i?S(s) ds, for 
i = 1,3, and G2 = PerQ - Girj. Equation (14.2 ip yields 

J(r) = ^Gir,r2 and $(r) = ^Ga P - 2ai$o(0] , (6-10) 

where <i>„(r) — J^f^fr^ds, for n = 0, 1. Therefore, for G3 = Gi?7 — G2, 

(J - $)(r) - ^ (Gsr^ + G2 2ai$o(r-)) . (6.11) 

Some algebra reveals that /,i(l) — {X/2Ki)rl [G3 In,i + G2 2ai In,2], where Ii^i — 

1 - ai, h,2 = $1(1) - 0.2, 13,1 = 2r^/i,i - ai, 13^2 = '"o [h.2 ~ ^h^i] + ^ - 02, and 

02 = ai$o(l). Hence, 

- + TJX ('^''^''»' + t5j:x?^'('»>) + (02) 

where, recalling that Tq = Ki/{2JiX/3ri), 
biirl) = 2ai(a2-l)+r^2ai(/ia-ai/i,2)], 

&2(r^) = ai(l - 2a2 + 2/3,2) - 2/3,1 + r2{ai[2/i,2(2ai - 1) - 3/i,i] + 2/i,i}, 
^3(^0) = /3,i-2ai/3,2+r2{ai[2/i,2(l-ai)+/i,i]-/i,i}. (6.13) 
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(ii) Drift and effective diffusivity, 77 ^ 1 

In this limit, the cells are affected by the flow to the extent that they mostly tum- 
ble. Using ([6?5]), equations (031), (|4?T9| and (|4J8| provides i?[J(r) = (1 - nj)r-^'^, 
where 137 = P/irj, mg(r) = r^C^-^) and 

A*o{r)=Per'''^ (r^-^^^^rA , (6.14) 



respectively. Hence, 



Ao=AS(l) = Pe-^. (6.15) 



Similarly, with definitions (|6.5p . equation (|4.2ip yields 



<I'M=Pe^(.^-y). (6.16) 
Hence, the equation (|6.2p for the effective diffusivity gives 

De = ^ + 2{Ao-Pc) [ r^R^dr + APef r^<i>R^dr = ^ +Pe^G{w), (6.17) 



where 



3 1_-w 
= — — 



1 — TU 

2 — 137X3 — t37 2 — UJ I 3 — n7 4 — 137 



(6.18) 



(iii) Dependence of dispersion on flow parameters in the strong and weak limits 

Here, drift and diffusivity are evaluated as a function of Pe for realistic param- 
eters. Recalling Pe = Ua/D", /3 = Vsa/D", 77 = UB/a and A = \/{2Bdr) (where 
dr is the rotational diffusion constant for swimming cells) we see it is in theory 
possible to vary Pe whilst holding /3, 77 and A (and so rg, zu) fixed. For C. nivalis 
the gyrotactic reorientation time B — 3.4 s, dr — 0.067 s^^, and so A = 2.2, thus 
Ki = 0.57, K2 = 0.16, Ji = 0.45, J2 = 0.16 (Pedley & Kessler 1990; HiU & Hader, 
1997). With these values, Gi = — 4(J2 — JiKi) = 0.39. Furthermore, the average 
swimming speed and cell diffusivity are Vg ~ 10~^ cm s~^ and D'^ 5 x lO""* cm^ 
s~^ (Hill & Hader 1997, Vladimirov et al. 2004). Using these parameters and a = 1 
cm, we find that /3 = 20 and rg = a/2 for 77 = 0.007, tq = 0.22 for 77 = 0.3, and 
m = /3/(477) = 0.05 for rj = 100. Hence, expressions dill]), ((02|) . ([QS]) and (l6T7)) 
are used to plot the effective diffusivity and drift (inset) for algae in a Poiseuille 
flow in figure 13^. The figure reveals that low and high levels of gyrotaxis (measured 
by 7/) lead to behaviour akin to Taylor dispersion but intermediate levels dramat- 
ically reduce the impact of advection. This is because at these intermediate levels 
the cells form dense plumes in the centre of the tube so are not subject to the full 
range of flow speeds. On the other hand, intermediate gyrotaxis does lead to large 
amounts of swimming and flow induced drift relative to the mean flow, due to their 
central location. For small and intermediate 77 the asymptotic results reveal that 
the drift changes sign for a non-zero Pe number. However, the asymptotic results 
for large 7;, are not strictly valid for small Pe. Nonetheless, one would expect the 
drift to change sign in a similar manner, such that all three curves intersect on the 
y-axis. 
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Figure 4. Effective diffusivity (inset: drift, Ao) against Pe calculated using 
asymptotic solutions to the Fokker-Planck approach, a) Poiseuille approxima- 
tion (case I) for r) = 0.007 (solid line), 77 = 0.3 (long dashed), r) = 100 (dashed) 
and the classical Taylor- Aris result (grey), b) Self-driven flow (p^ = 0) using 
the coupled solution from case II (broad plumes; rj — 0.007). Mode 1 (long 
dashed) and mode 2 (dashed) are shown with the uncoupled limit (solid). 

(d) Algae in self-driven flow (weak coupling, A<S^ 1) 

Recall from section [3] that for self-driven flows there are two solutions: a simple 
mode 1 flow and a mode 2 flow with upwelling at the tube sides. To calculate the 
transport coefficients in these cases the definitions (|6.4p are employed (since the flow 
solutions were all obtained for 77 ^ 1 and weak coupling, A <^ 1). With the same 
parameters as for the 77 = 0.007 case in figure H^, figure HId plots the diffusivity and 
drift (inset). It is clear that the mode 1 results for these broad plumes are rather 
similar to those generated by Poiseuille flow. However, for the mode 2 solutions 
cells both drift and diffuse faster, likely due to the greater shear. 

7. Discussion 

In this paper, we derive exact expressions in the long-time limit for the mean 
drift and effective axial diffusion of an axisymmetric blob of biased, swimming 
microorganisms in a plume in a pipe flow driven by an external pressure gradient 
and the presence of the (negatively) buoyant cells. In the same limit, we find that 
the axial skewness of the cross-sectionally averaged cell distribution vanishes. The 
results are independent of the cell geometry, swimming behaviour and model used 
to represent the cell-flow interactions. 

Explicit results for several useful cases are presented, from the Taylor- Aris limit 
to fully coupled gyrotactic spherical swimming cells (i.e. cells that drive the flow 
and whose swimming direction is biased by external and viscous torques). The ex- 
pressions reveal the mechanisms for several competing effects and explain how these 
lead to diffusion and (positive or negative) drift through the tube. Fundamentally, 
the cells swim and, in the limit that they are very bottom heavy, they may swim 
mostly against a downwelling flow, leading to a negative drift relative to the mean 
flow. On the other hand, cells that are not bottom heavy act more like diffusing 
passive tracers, with no drift. In both these cases the cells diffuse as for Taylor-Aris 
dispersion. However, an intermediate degree of bottom heaviness leads to much 
more interesting behaviour. A balance between gravitational and viscous torques, 
a balance that will vary across the pipe flow, can lead the cells to form gyrotac- 
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tic plumes, inducing further flow and self-concentration. These centrally focTiscd 
plumes of cells can be strongly advected with the flow (i.e. faster than the mean 
flow) but will sidestep classical shear-induced Taylor- Aris dispersion; effective dif- 
fusion may be dominated by swimming diffusion even for large flow rates. It is clear 
that swimming behaviour leading to drift across streamlines can have a tremendous 
influence on cell transport in such systems. 

The results are sufficiently general that they may easily be applied to other 
micro-organisms and taxes, such as chcmotaxis in suspensions of bacteria swimming 
in flows in microfluidic chambers, or spermatozoa in vivo. In a subsequent paper we 
shall provide further explicit examples for non-spherical cells (behaviour influenced 
by the rate-of-strain tensor) and for additional swimming stresses for concentrated 
suspensions. Both these aspects will modify the plume structure and thus affect 
axial cell transport. 

Work in progress is exploring how the theory can be applied to determine 
the qualitative form of the orientationally averaged cell swimming diffusion tensor 
for suspensions of gyrotactic cells from experiments. For a realizable experiment 
one must introchicc dyed cells into a plume whilst maintaining a constant cross- 
sectionally averaged cell concentration. This may be achieved simply by momen- 
tarily switching from undyed to dyed cells at the input or using photoactivatable 
GFP for localized photolabelling of cells (Patterson & Lippincott-Schwartz 2002). 
Note that plume solutions for the various diffusion descriptions differ qualitatively 
for large Peclet numbers, and thus so must predictions for mean drift and effective 
diffusion. Hence, we aim to clarify the applicability of differing diffusion approxi- 
mations in a general shear flow. 
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